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1. Introduction 

The theory of electromagnetic wave propagation in dispersive nonhnear media such as 
optical fibers is one of the key elements of the information age in the 20th century. 
Maxwell's equations present a complete description of the classic (non-quantized) 
electromagnetic field. In the presence of nonlinear effects, however, a general solution to 
Maxwell's equations seems impossible and approximate solutions are commonly used to 
explain phenomena. For the particular case of pulse propagation in cubic nonlinear 
media, the cubic nonlinear Schrodinger equation (NLSE) presents a slowly-varying 
amplitude approximation that has been extremely successful in nonlinear optics [T]. 
The NLSE belongs to the small class of integrable partial differential equations that can 
be solved by inverse scattering transform. It admits solitary wave solutions that show 
stable propagation in Maxwell's equations as long as one is working in regimes where 
solutions of the NLSE are close to solutions of Maxwell's equations. The extraordinary 
stability properties of optical solitons over long distances have made them promising 
candidates for bit carriers in fiber-optic communications jTO] . 

The derivation of the NLSE assumes a small parameter which is the ratio of pulse 
width and period of the carrier wave. This small parameter is used in creating an 
asymptotic expansion, typically using multi-scale expansions or the renormalization 
group method. In the regime of ultra-short pulses, where the pulse lengths shorten, 
the NLSE approximation becomes less accurate [HI [2l]. It is possible to increase 
the accuracy of the NLSE approximation by incorporating higher-order corrections. 
Physically, these next-order correction terms model important effects as Raman delay, 
higher-order dispersion, or self-steepening of the pulse pQ. 

Advances in the generation of ultra-short pulses allow pulse widths in the range 
of 10 fs or below. The nonlinear optical pulse compression technique even makes it 
possible to generate pulses of about 6 fs pT]. Femtosecond laser systems generating 
short pulses have led to diverse interests, studies and apphcations of ultra-short pulses 
such as extremely high intensity laser-matter interactions, characterization of high-speed 
electronic and optoelectronic devices and systems, optical communications, medical 
imaging as well as ultrafast phenomena in solid-state, chemical and biological materials 
[25]. 

Recently, the so-called short-pulse equation (SPE) was proposed [23] to describe 
nonlinear pulse dynamics exploiting a scaling designed to describe the behavior of pulses 
in the ultra-short limit. Sakovich and Sakovich showed that the SPE is integrable [19] . 
and possesses exact one-solitary wave [20], later multi-solitary wave solutions [T3] were 
found as well. This development led to an intensive research on the SPE over the years. 
The bi-Hamiltonian structure [2], [3] , conserved quantities [9] and the periodic solutions 
[I6] of the SPE were analyzed. Generalizations as the vector short pulse equation [T7[l22] 
and the regularized short pulse equation (RSPE) [8J were studied. For the RSPE, the 
existence of multi-pulses [2], traveling waves, and solitary wave solutions [SI |7] were 
confirmed. 
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It is natural to ask whether the sohtary wave solutions of SPE persist in Maxwell's 
equations when chosen as initial conditions in the ultra-short pulse regime, similar to 
the sohton solutions of NLSE in the broad pulse regime. Numerical studies show that 
this is indeed the case [13]. Observed over longer distances, however, the solitary wave 
solutions of SPE show slow distortions when propagating in Maxwell's equations. This 
behavior is expected, as the SPE is only an approximation to Maxwell's equations. 
Similar to the NLSE, it is therefore worthwhile to consider higher-order contributions 
in order to obtain a better approximation to Maxwell's equation. 

The derivation of such higher-order corrections to the SPE is the objective of 
the present paper. In the next section, following the previous work we use the 
renormalization group method in order to derive the higher-order short-pulse equation 
(HSPE). Section [3] presents a detailed numerical study of the new equation, comparing 
results to the original SPE and to Maxwell's nonlinear wave equations. 



2. Derivation of the higher-order SPE 

In this work, we analyze a one-dimensional nonlinear wave equation describing the 
propagation of pulses in nonlinear cubic media given by 

Uxx = utt + XoM + xA{'^^)tt\ (1) 

where the Xo and Xs model the response of the medium to the applied field. Eq. ([T]) is 
derived directly from Maxwell's equations under certain assumptions, for details we refer 
to previous work [13]. In this work, we assume the operators xo and X3 of particularly 
simple form, namely as factors in Fourier domain, meaning that if u denotes the Fourier 
transform of m, 

Xk[u] = XkU ■ H {\uj\ - ^/Xa) ^ fc = l,3. (2) 

(with "H being the Heaviside-function) such that unstable linear modes of ([T]) are 
excluded. The numbers Xk characterize the strength of the response of the medium. 
We are interested in asymptotic solutions of the form 

ij — X 

M(x,)f:) = eAo(0, Xi,X2, ...) + e^Ai + 0= — - — , a;„ = e"x. (3) 

and, it has been shown [23] that, at the leading order, the evolution equation for Aq is 
given by 

- 2d^dx, = Xo^o + XzdlAl . (4) 
or, more precisely as 

- 2d^dx,A^ = UM + X3 [dlAl\ (5) 

where the cut-off in the operators xo and xs has been re-scaled appropriately, hence the 
cut-off frequency being now ey/xo- the following, for simplicity of notation, we will 
write the susceptibilities as factors instead of operators. 



Higher-order corrections to the short-pulse equation 



4 



We also note that equation Q can be written after a simple transform as [2T] . 

UxT = U + Vix • (6) 

In the following, we compute higher-order corrections of this expansion. Although it is 
possible to do this using a multi-scale expansion \l3l , we prefer to present the derivation 
using the renormalization group (RG) method. The RG method, which was developed 
as a tool for asymptotic analysis [H |5], has already been successfully applied to the 
nonlinear wave equation ([T]) to derive the SPE Q [6] at the leading order. This method 
has two main advantages over the method of multiple scales: (a) algebraic calculations 
in the case of the RG treatment are simpler than the ones appearing in the multiple-scale 
expansions especially if higher-order corrections are considered, and (b) the introduction 
of a particular multi-scale- ansate is not required a priori for the reason that the RG 
equation naturally gives rise to the scales that need to be introduced for a consistent 
asymptotic expansion. 

We first reproduce quickly the SPE and then proceed to the next order. We apply 
a coordinate transformation in the form of 

u{x,t) = B{(j),x) (7) 

where and x are given in ([s]). Substituting ([t]) in ([T]), we obtain 

2 1 
— -B^x = XoB — Bxx + -^XsiB . (8) 

We assume the expansion 

5(0, x) = eAo(0, x) + e^Aiicp, x) + ■ ■ ■ . (9) 

Substituting the expansion ^ into ([s]), the terms of 0(1) yield 

-2(Ao)<^. = 0. (10) 

This implies Aq is independent of x, i.e., Aq = Ao(0). On the other hand, the terms of 
0(e) yield 

-2(Ai)<^, = XoAo + X3(Ao)<A</.- (11) 
Therefore, we find the second order approximate solution 

B^'\<P,x) = eAo(0) - ( r XoAo(0)# + X3 (Ao(0)')^) • (12) 



Notice that a secular term appears on the equation (12) which corresponds to the 
term proportional to x. In order to remove this secular term, we consider the term 
Ao(0) — f a^(/^oo Xo Ao(0)c^0+X3 (Ao(0)^)<)<,) as the Taylor expansion of a function Vo(0, x) 
about X = 0. Thus, we need to find Vo(0, x) which satisfies 

Vo(0,x = O)=Ao, (13) 



dx 2 



XoVo(0)# + X3 (Vo(0)')j . (14) 
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This equation introduces a new scale ex. Let us define xi = ex as in equation ([s]) then 
we obtain 



Vo(0,xi =0) =Ao 
dxi 2 



XoVo(0)# + X3 (Vo(0)')^ 



(15) 
(16) 



Solving equation (16) provided that the initial condition ([15j) is satisfied, we can express 
the r.h.s. of equation (12) as the Taylor expansion of Vq. Hence, we finally obtain the 
second order approximate solution, B^'^\(f),x) 

1 



("1^0)0x1 



2X0V0 - - 



eVo- Furthermore, equation (16) yields 

(17) 



This is the SPE. 

For higher order corrections, we follow the similar steps. First, we need to collect 
higher order terms in e. The usual way of obtaining these is by assuming the ansatz 
B = eAo + e^Ai + e^A2 + - ■ ■ and collecting appropriated terms. However, this will lead to 
highly complicated algebraic calculations. Here, we approach this problem by assuming 
a different ansatz. Since we have already obtained the second order approximation of 
the solution, B^'^\ we assume 

5(0,x) = 5(2) + e'=+iAfc(0,x). (18) 



Plugging the equation (18) into equation ^8^, we find that the next higher-order 



nontrivial contribution comes from the term —B^x in (jSj), corresponding to — e^(Vo)a;ixi- 
Hence we need to choose A; = 3. Then, O(e^) terms yield 

- 2(A3)^. = -(Vo).,..i . (19) 

Before we proceed with (19), let us define A^(Vo) = — (Vo)xia:i and rewrite this term. 
From the SPE (|l7|), it follows that 

1 



-(V, 



Ojxi 



Vom~^+\xz (Vo(0)')^. 



Therefore, using the second derivative of Vq, 

- (Vo).,., = \xof_ {Vo)xAm+\x3 {^o{ct^f)^x, ' 
and the SPE along with some algebraic manipulation, we obtain A^(Vo) 



(20) 



(21) 



A^(Vo) 



4 



3 



# / Vo(0')ci0' 



:XoX3Vo(Vo)<^ / Vo(0)rf0 



- 2X3Vo(Vo)^(Vo^; 



X3Vo'(Vo^), 



(22) 
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Clearly, A^(Vo) involves higher-order linear and nonlinear terms, a more detailed 
discussion will follow in the next section. Coming back to (19), we obtain the 4th 
order approximate solution 



B^^\(P,x) = e Vo(0,xi) - -X / iV(Vo)# 



(23) 



In order to combine all the possible secular terms, we rewrite equation (23) using 
equations (13),(14) and Taylor expansion of Vq. This yields 

i?(^)(0,x) = e fAo(0) -^x( xoAo(0)# + X3 

(24) 



iV(Ao)#' + 0(e 



Again, to remove the secular terms which are proportional to x, we now need to find 
A{(j), x) satisfying 

A(0,a; = O)=Ao(0), (25) 



dA 

dx 



XoA(0)rf0 + X3 (^(0)'), 



N{A)d(t^ . (26) 



XoAi4>)d^ + X3 (^(0)'), 



N{A)d(j)' . (27) 



The term that is proportional to e'^ leads us to introduce a new scale xs = e^x in addition 
to Xi = ex. From equation (26), we also find 

dA 

dx~i~'2 
Denoting x = xi, we finally obtain 

- 2A^^ = xoA + X3(^')00 + e^N{A) . (28) 

This is the higher-order short pulse equation, and A{>c, (p) is the magnitude of the electric 
field following the introduction of the new variable x. It is obvious that this equation is 
an extension of the SPE ^ with additional terms on the right hand side. The additional 
operator N{A), represents the higher-order corrections of the SPE approximation. 



3. Numerical Analysis 

In order to show that the HSPE captures the effects beyond the SPE approximation, 
we compare the solutions of the original nonlinear wave equation ([T]) to the solutions 
of the SPE and the HSPE. As initial conditions, we use one- and two-soliton solutions 
of the SPE. We first summarize very briefiy the structure of these solutions and then 



discuss their evolutions in the HSPE given by (43) in comparison to the nonlinear wave 
equation ([T|. Both equations are numerically solved using standard techniques. For the 
nonlinear wave equation we use a pseudo-spectral leap-frog scheme in Fourier space, the 
SPE and the HSPE are solved using a fourth-order exponential time-differencing method 
[12] in Fourier space. For the solution of the nonlinear wave equation it is essential to 



implement the cut-off in the susceptibility operators as in (g. For the HSPE, this can 
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be done in a similar fashion. In the examples considered here, we have found only a 
very little difference when neglecting the cut-off and treating the operators simply as 
scalars. However, for consistency, the cut-off was implemented in all simulations. 

Before we proceed, let us also note that, for the numerical simulation, it is useful 
to cast the HSPE in a slightly different form, 

A. = -f l'^A{>.A)d<P-f {A')^-'^ l\{A{Kj))d4>, (29) 

where we write the integral of the operator N{A) involving the higher-order correction 
terms as 

AO 



N{A{x,(P))d(P=^ d<Pi d(f)2 d(j)3A{>c 

^ J ~oo J —oo J —oo 



4 

3xoX3 
4 

-f^A(x,0f (A(x,0f), . (30) 



A(x, 0)^ / (i0A(x, 



,2 f4> f4>i f4>2 

AO 



Here, we define 

L,{A) = ^ I d<f)i I d<j)2 I d<f)sA{>c,<j)3), (31) 

^ J~oo J — oo J —oo 

N,iA) = ^ r dMi^,4>?, (32) 
N,(A) = ^A(x, 0)2 r dMi^, 0), (33) 

^ J — oo 

iV3(A) = ^A(x,0)2(A(x,0)=^)^. (34) 

In the next subsection, we will compare the linear term Li to the nonlinear terms 
Ni, N2, N3. Also note that numerical evaluation of these terms can be done by using 
the Fast Fourier Transform and that the singularity at a; = is removed by introducing 
the cut-off. 



3.1. Propagation of one- pulse solutions 

The analytical one-soliton solution of the SPE is found 

m sin -0 sinh + n cos cosh 

^ = ^"'^ 2 ■ 2 , ^ 2 IJ^ (2^) 

m"^ sm -0 + 71^ cosh 



with 



m sin 2ib — n sinh 2(w , , 

X = Y + 2mn ^ . / (36) 
sm -0 + cosh 



m(F + r), i) = n{Y-T), n= ^{l-m?). (37) 
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Here, we only consider non-singular solutions, which implies that the soliton parameter 
m satisfies the condition < m < sin | ^ 0.383. 

In our numerical simulations, we set m = 0.35 in order to have strong nonlinear 
effects, but no singularities. We also choose e = 0.4 since we are mainly interested in 
higher-order effects. Therefore, we expect that the solutions of the nonlinear wave 
equation and the SPE will show differences already after fairly short propagation 
distances. 




Figure 1. Snapshot of the evolution of the SPE soliton (35) in the nonlinear wave 
equation ^ (soHd hne), the HSPE Q (dahed hne) and the SPE ^ (dotted hne). 



We set a propagation distance x = 51.2, and susceptibilities xo = 2, xa = 1/3. 



Using an SPE soliton (35) as the initial condition, we numerically integrate the SPE 



(|6|, the HSPE (|43|, and the full nonlinear wave equation (ij). Figure [l] demonstrates the 
numerical solutions at the distance x = 51.2 of the SPE, the HSPE, and the nonlinear 



wave equation. Clearly, if an SPE soliton (35) is taken as the initial condition, the 



HSPE provides a much better approximation to the nonlinear wave equation than the 
SPE. 



We now consider the nonlinear operator (22) that presents the higher-order 
contributions. Although, we cannot make a general statement, our numerical 
experiments show that for large m-values close to singularity, surprisingly, the linear 
term seems to be dominant. 

Figure [2] shows all four terms, Li, Ni, N2, and N3 for the Sakovich-soliton with 
m = 0.35 and, clearly, the nonlocal dispersive term is much stronger than the nonlinear 
terms. This numerical observation leads to a truncated higher-order SPE, that only 
involves the linear term and, therefore, is much easier to compute and to analyze: 



2 Xo 



dcj) / A(0/)c/0 



(3J 
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Figure 2. Comparison of the four correction terms from eq. (301, clockwise starting 
at the upper left figure: Li, Ni, N2, and N3, for a Sakovich soliton as initial condition 
at a; = 0. 



Figure [3] shows the numerical solutions of the HSPE and the truncated HSPE in 
comparison to the numerical solution of the SPE, again with a Sakovich soliton being the 
initial condition. This figure shows that the difference in the performance of truncated 
HSPE and HSPE is very small and the truncated HSPE is still a better approximation to 
the full nonlinear wave equation than the SPE. This implies that, for practical purposes 
in some cases, the truncated HSPE might provide sufficient additional accuracy. 

3.2. Multi-soliton propagation 

Multi-soliton solutions of the SPE have been derived [IS] through a systematic procedure 
from breather solutions of the sine-Gordon equation. The parametric multi-soliton 
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3 




-0.2 



Figure 3. Snapshot of the evolution of the Sakovich soliton in the HSPE (solid line), 
the truncated HSPE (dashed Une) and the SPE (dotted line). 



solution can be expressed in the compact form as 



with 



uiy,t) = 2i(^ln^^ , xiy,t)=y-2(lnf'f^^ + d 



^=0,1 j=l l<j<N 

N 



M=0,1 3=1 

= PjV + —t + CjO, 

Pj 



l<j<N 

(j = l,2,...,M) 
{j,k = l,2,...,M;j^k) 



(39) 



(40) 



Pj + Pk 

P2j~i = P*2j = aj + ibj, aj > 0, bj > 0, {j = 1, 2, M) 

6j-i,o = C2j,o = + (i = 1, 2, M) 

9j = aj{y + Cjt) + A,-, {j = 1, 2, M) 
Xj = - Cjt) + Hj, {j = 1, 2, M) 

^^• = ;;2^' (j = i,2,...,M), 

where pj and are arbitrary parameters such that pj ^ ±pk ior j k, i — y/—!, N 
is an arbitrary positive integer, and M — N/2 is the number representing the multi- 
sohton solutions (one soliton, two solitons, etc.). If = 4 and M = 2 are chosen, 
one can generate the two-soliton solution. The condition for a single- valued nonsingular 
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multi-breather solution is 

M 



(41) 




Figure 4. Snapshot of the evolution of the SPE two-soUton solution (39) in the 



nonlinear wave equation ([T]) (sohd hne), the HSPE (43) (dashed line), and the SPE 
(dotted line). 

In order to investigate collisions of SPE-solitons in the HSPE and the nonlinear 

wave equation, we use a two-soliton solution as the initial condition. The parameters 

are chosen ai = 0.1, hi = 0.5, a2 = 0.16, 62 = 0.8, d = 0, Xi = 10, A2 = 0, /ii = 0, 

/i2 = 0. Here, the propagation distance needs to be much longer in order to capture 

the entire collision. Therefore, we choose e = 0.2 and a propagation distance x = 375. 

At such large distances, higher-order terms start to influence the propagation as it can 

be seen in figure |4} Clearly, the HSPE approximates the solution of the nonlinear wave 

equation much better than the SPE. However, it seems to be necessary to incorporate 

even higher order terms in the approximation derived in section [2j Extending the results 

of section |2| it is easy to see that the next linear term in the expansion is given by an 

integral operator involving four integrations of the form 

^3 /•<j) f4>i r4>2 r4>i 

e^M{A)=e^f d(f)J d(t>2 d(f), d(j),A{x,(j),) (42) 
^ J —00 J —00 J —00 J —00 

yielding the equation 

- 2A.^^ = xoA + xM'')h + e'iV(^) + e^M{A) . (43) 

In a heuristic way, this additional term can be understood in the following way: 
Consider again equation ([s]), we see that higher-order contributions will arise from —B^x 
in dsl), in particular a term —2Ax^x:i which will be renormalized by introducing a scale 
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X5. From equation (43), we find that 2Ax.i^ 
terms, we obtain 



A 



XlXl • 



Ax. 



02 



,A 



Therefore, keeping only hnear 

(44) 



and hence, 



A ~ 



16 



02 



iA 



(45) 



yielding the addtional linear term in the HSPE (43). 

Indeed, figure |5] shows that there is a considerable improvement for the evolution of 
two-soliton solution when incorporating this additional linear term on the r.h.s. of eq. 



(43). Without showing the corresponding figure, we also note that, as expected, there 



is improvement for the approximation of the one-soliton solution. 
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Figure 5. Snapshot of the evokition of the SPE two-sohton sohition (391 in the 
nonUnear wave equation ([T]) (sohd Une) , the HSPE ( 43 1 with the additional hnear 
term (42) (dashed hne) and the SPE (dotted hne). 



4. Conclusion 

Using the renormalization group method, we derived higher-order correction terms to the 
short-pulse equation. We showed numerically that the incorporation of these terms can 
considerably improve the accuracy of the SPE approximation when describing solutions 
to Maxwell's equations. This was shown using one- and two-sohton SPE solutions as 
initial conditions. Even the incorporation of one additional linear term in the SPE 
approximation can lead to considerable improvement, and the extension to higher-order 
terms is straight-forward. A precise study of the effect of the higher-order terms on the 
solitary waves is subject to future studies. 
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